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Abstract 


In  the  context  of  the  surveillance  of  space,  known  resident  space  objects  (RSO),  i.e.  active 
satellites  and  space  debris,  need  to  be  acquired  periodically  to  ensure  that  the  knowledge  of  their 
orbital  parameters  is  up-to-date.  The  most  sensitive  observation  method  consists  of  acquiring  an 
image  sequence  while  the  target  is  actively  tracked  by  the  optical  system  such  the  RSO  appears 
stationary  in  the  image  sequence.  This  generates  an  incredible  amount  of  data  which  requires 
efficient  detection  algorithms  and  automated  processing  software.  In  this  context,  this  document 
presents  models  for  image  formation,  acquisition,  and  processing.  From  the  observation  request 
to  the  production  of  the  images,  the  stream  of  events  is  briefly  modeled.  This  ensures  that  a-priori 
knowledge  and  strategy  will  be  applied  to  efficiently  process  the  acquired  data.  These  models 
produce  a  series  of  algorithms  used  to  separate  the  image  components  and  detect  the  RSO.  This 
document  presents  algorithm  diagrams  that  can  be  used  to  develop  the  processing  software. 


Resume 


Dans  le  cadre  de  la  surveillance  de  Tespace,  les  objets  spatiaux  connus  en  orbite  (OSO),  i.e., 
satellites  actifs  ou  debris  spatiaux,  ont  besoin  d’etre  acquis  periodiquement  pour  s’ assurer  que  la 
connaissance  de  leurs  parametres  orbitaux  soit  maintenu  a  jour.  La  methode  d’observation  la 
plus  sensible  consiste  a  acquerir  une  sequence  d’images  lorsque  la  cible  est  poursuivie  activement 
par  le  systeme  optique  et  que  l’OSO  apparait  stationnaire  dans  la  sequence  d’images.  Cela  genere 
une  quantite  incroyable  de  donnees  qui  requierent  des  algorithmes  de  detection  efficaces  et  des 
logiciels  de  traitement  automatise.  Dans  ce  contexte,  ce  document  presente  des  modeles  pour  la 
formation,  Tacquisition  et  le  traitement  des  images.  De  la  requete  d’observation  a  la  production 
des  images,  la  chaine  des  evenements  est  modelisee  brievement.  Cela  permet  de  voir  les 
connaissances  et  strategies  a  priori  qui  peuvent  etre  utilisees  pour  traiter  efficacement  le  donnees 
acquises.  Ces  modeles  produisent  une  serie  d’algorithmes  utilises  pour  separer  les  composantes 
de  l’image  et  detecter  l’OSO.  Ce  document  presente  les  diagrammes  d’algorithmes  qui  seront 
utilises  pour  developper  le  logiciel  de  traitement. 
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Executive  summary 


Image  and  processing  models  for  satellite  detection  in  images 
acquired  by  Surveillance-of-Space  Space-based  sensors 

Martin  P.  Levesque;  DRDC  Valcartier  TR  2009-095; 

Defence  R&D  Canada  -  Valcartier;  . 

In  Surveillance  of  Space  (SofS),  resident  space  objects  (RSO)  must  be  periodically  observed  and 
the  measurement  of  their  position  is  used  to  update  their  orbital  parameters.  This  huge  task  is 
assumed  by  the  Space  Surveillance  Network  (SSN),  a  component  of  NORAD.  Canada 
contributes  to  these  measurements  with  several  small  telescopes.  Two  DND  satellite  projects  are 
also  in  development  for  this  purpose.  They  are  the  SAPPHIRE  and  NEOSSat  satellites. 

This  document  contains  models  that  describe  the  image  formation  and  acquisition  process  for 
sensors,  ground-based  or  space-based,  that  acquire  sequences  of  images  of  these  RSO.  Based  on 
these  models,  it  describes  processing  strategies,  algorithms  and  diagrams  for  the  automatic 
detection  and  reporting  of  RSOs.  This  document  succinctly  presents  the  lessons  learned  from  the 
preceding  multi-year  joint  (DSpaceD-DRDC)  effort  on  satellite  detection  using  images  acquired 
in  sidereal  tracking  mode.  This  study  applies  these  lessons  learned  to  detection  using  image 
sequences  acquired  in  track  rate  mode  (TRM).  TRM  is  the  most  sensitive  method  and  is  the 
acquisition  mode  that  will  be  implemented  in  the  SofS  satellites  in  development. 

This  study  is  a  contribution  to  the  data  exploitation  of  the  NEOSSat  satellite.  The  development  of 
models  and  algorithms  is  the  starting  point  for  the  conception  of  data  processing  software.  This 
document  provides  guidelines  for  software  development.  Even  if  the  SAPPHIRE  exploitation 
software  is  already  in  development,  this  document  contains  guidance  on  algorithms  and 
processing  strategies  that  can  improve  detection  performance  considerably.  The  proposed 
processing  method  may  be  used  to  upgrade  the  processing  software  for  this  satellite. 

This  processing  model  is  designed  for  post-acquisition  processing  by  ground  stations  where  all 
the  required  computing  resources  are  available.  It  is  not  expected  the  required  processing 
capability  is  available  on  board  the  satellite. 
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Dans  le  cadre  de  la  surveillance  de  l’espace  (SdeE),  les  objets  spatiaux  en  orbite  (OSO)  doivent 
etre  observes  periodiquement  et  la  mesure  de  leur  position  est  utilisees  pour  mettre  a  jour  leurs 
parametres  orbitaux.  Cette  tache  enorme  est  assumee  par  le  ‘Space  Surveillance  Network’  (SSN), 
une  composante  de  NORAD.  Le  Canada  contribue  a  ces  mesures  grace  a  plusieurs  petits 
telescopes.  Deux  projets  du  MDN  sont  aussi  en  developpement  pour  cette  tache.  Ce  sont  les 
satellites  SAPPHIRE  et  NEOSSat. 

Ce  document  contient  des  modeles  qui  decrivent  la  formation  des  images  et  le  processus 
d’acquisition  de  capteurs,  bases  au  sol  ou  dans  l’espace,  qui  acquierent  des  sequences  d’images 
de  ces  OSO.  En  se  basant  sur  ces  modeles,  il  decrit  des  strategies  de  traitement,  des  algorithmes 
et  des  diagrammes  pour  la  detection  et  la  declaration  automatique  des  OSO.  Ce  document 
presente  succinctement  les  Iccons  apprises  de  plusieurs  annees  d’ efforts  communs  (DspaceD- 
DRDC)  sur  la  detection  de  satellites  utilisant  des  images  acquises  en  mode  de  poursuite  siderale. 
Cette  etude  applique  ces  Iccons  apprises  a  la  detection  en  utilisant  des  sequences  d’images 
acquises  dans  le  mode  de  poursuite  des  satellites  (TRM  :  track  rate  mode).  C’est  la  methode  la 
plus  sensible  et  c’est  le  mode  d’acquisition  qui  sera  implante  dans  les  satellites  en  developpement. 

Cette  etude  est  une  contribution  a  l’exploitation  des  donnees  du  satellite  NEOSSat.  Le 
developpement  des  modeles  et  des  algorithmes  est  le  point  de  depart  pour  la  conception  du 
logiciel  de  traitement  des  donnees.  Ce  document  sera  utilise  comme  guide  pour  ce 
developpement.  Meme  si  le  logiciel  d’exploitation  de  SAPPHIRE  est  deja  en  developpement,  ce 
document  contient  des  algorithmes  et  des  strategies  de  traitement  qui  peuvent  ameliorer 
considerablement  la  performance  de  detection.  La  methode  de  traitement  propose  peut  etre 
utilisee  pour  une  prochaine  mise  a  niveau  du  logiciel  de  ce  satellite. 

Ce  modele  de  traitement  est  concu  pour  le  traitement  differe  par  les  stations  au  sol  ou  toutes  des 
ressources  informatiques  sont  disponibles.  On  ne  s’attend  pas  a  ce  que  la  capacite  de  traitement 
requise  soit  disponible  a  bord  du  satellite. 
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Introduction 


To  maintain  a  good  level  of  space  situation  awareness,  the  orbital  parameters  of  the  Resident 
Space  Objects  (RSOs;  satellites  and  space  debris)  require  frequent  updates.  Most  of  the  RSOs 
gradually  deviate  from  their  known  orbits.  The  main  causes  of  the  deviations  are  the 
unpredictable  solar  activities  with  the  resulting  variable  atmospheric  drag,  and  the  irregularities  in 
the  Earth’s  gravitational  field.  Thus,  the  RSOs  require  periodic  observations  and  accurate 
position  measurements  for  the  recalculation  of  their  orbital  parameters.  Most  of  these 
observations  are  done  with  ground-based  telescopes,  but  space-based  observation  systems  provide 
more  consistent  performance  (no  meteorological  limitations,  potential  for  24  hour  per  day 
operation)  and  are  becoming  more  available. 

Because  of  the  quantity  of  RSOs,  all  the  update  processes  (tasking,  observation  and  data 
processing)  must  be  automated  as  much  as  possible.  Usually,  the  degradation  rate  of  specific 
orbits  is  known  and  the  tasking  process  can  be  automated  by  selecting  the  TLEs  older  than 
predetermined  values.  A  TLE,  which  means  ‘Two  Lines  Element  sef,  is  a  database  table  format 
that  contains  the  orbital  parameters.  Then  a  scheduling  program  transforms  these  task  requests 
into  a  series  of  observation  requests  where  the  selected  sensors  should  have  the  optimal 
conditions  (time,  distance  of  RSO,  phase  angle,  etc.)  to  perform  their  observations.  The  sensors 
are  programmed  and  automatically  perform  their  observations,  acquiring  an  incredible  amount  of 
data.  These  data  are  then  processed  with  satellite  detection  algorithms  and  the  detections  are 
reported.  Finally,  the  TLEs  are  recalculated  and  the  database  updated  based  on  detection  data. 
This  overall  process,  which  required  considerable  human  effort,  can  now  be  almost  completely 
automated.  It  requires  only  minimal  human  intervention  for  system  maintenance  and  system  data 
integrity  verification  and  validation. 

This  document  addresses  the  automatic  processing  issue.  In  previous  work  (Refs  1  and  2), 
algorithms  were  developed  for  the  detection  of  satellite  streaks  in  astronomical  images.  These 
algorithms  are  adapted  for  images  acquired  by  ground  telescopes  operating  in  sidereal  tracking 
mode.  This  document  considers  the  processing  of  images  acquired:  1)  in  satellite  tracking  mode; 
and  2)  by  a  moving  space-based  sensor.  It  describes  an  image  acquisition  model,  which  leads  to 
the  following  processing  approach  and  supports  its  development. 

The  goal  of  this  document  is  to  provide  the  physical  and  mathematical  models  and  guidelines 
required  for  the  development  of  future  processing  software  that  will  be  used  with  the  coming 
space-based  surveillance-of-space  systems.  The  software  is  not  yet  developed  but  most  of  the 
algorithms  developed  for  the  ground-based  telescope  are  reusable.  So,  there  is  no  risky  black  box 
in  the  following  mathematical  demonstration.  The  fuzzy  function  boxes  like  the  ‘streak  detector’ 
or  ‘star  identification’  make  reference  to  available  (and  already  tested)  commercial  software  or 
algorithms  already  developed  and  tested  in  previous  works. 

The  work  presented  in  this  document  begins  with  the  development  of  acquisition  and  image 
formation  models.  This  allows  seeing  the  a-priori  knowledge  that  can  be  used  by  the  processing 
algorithms.  Then,  the  image  model  is  formulated,  indicating  the  image  components  that  can  be 
separated  out.  The  following  algorithm  development  systematically  extracts  each  one  of  these 
components,  until  it  is  able  to  isolate  and  report  the  RSO  signal.  Finally,  the  processing  scheme 
indicates  how  the  algorithms  are  linked  together  along  with  their  dependencies,  which  provides  a 
clear  view  of  the  global  processing  methodology. 
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2  Tasking  and  acquisition 


Once  a  target  has  been  designated  for  an  acquisition  task,  programs  are  executed  to  determine 
when  and  where  the  observation  conditions  (mostly  based  on  the  brightest  conditions)  will  be 
best.  The  ‘Two  Lines  Element  set’  (TEE),  which  contain  the  orbital  parameters,  are  the  input 
data  for  a  program  called  ‘the  propagator’,  which  calculates  orbits  and  predicts  the  satellite 
position  ‘Rp’  and  motion  at  the  required  time  of  observation.  This  operation  produces  a  table 
where  the  sidereal  coordinates  ‘(Ra,  Dec,  t)’,  i.e.,  right  ascension,  declination,  are  listed  as  a 
function  of  the  time.  Note  that  this  table  is  generated  for  specific  observer’s  positions  which  are 
known  at  every  moment  by  the  propagator;  this  observer  may  be  another  moving  satellite.  The 
‘propagator’  function  is  usually  assumed  by  COTS  (commercial  off  the  shelf)  software  like 
‘Satellite  Tool  Kit’  or  ‘TheSky’. 

At  the  moment  of  the  acquisition,  these  (Ra,  Dec,  t)  coordinates  are  used  to  generate  system 
control  commands.  At  the  appropriate  moment,  the  optical  system  is  pointed  at  the  expected  RSO 
position  and  the  tracking  rate  is  adjusted  to  maintain  it  in  the  FOV.  At  the  appropriate  time,  the 
camera  acquires  the  images.  Because  the  pointing  system  never  points  exactly  at  the  designated 
position  (but  with  an  acceptable  error),  the  measured  pointing  parameters  (Ra,  Dec,  Rotation 
and  t)  are  recorded  in  the  pointing  data  file  or  in  image  header. 

Hence,  at  the  end  of  the  acquisition  process,  the  available  data  are:  the  image  sequence,  the  table 
of  expected  satellite  positions  and  the  measured  pointing  parameters.  In  addition,  the  image 
processing  scheme  can  use  a  star  catalogue  for  telemetry  evaluation  (astrometry)  and  the 
knowledge  of  the  acquisition  system  (field  of  view,  PSF  width,  etc.).  This  tasking  and  acquisition 
process  is  illustrated  in  Fig.  1. 
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Figure  1  -  Tasking  and  acquisition  sequence  showing  the  a  priori  knowledge  and  the  generated 

data 
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3  Scenario  of  image  acquisition 


The  RSO  can  be  acquired  in  two  different  modes;  Sidereal  Stare  Mode  (SSM)  or  Track  Rate 
Mode  (TRM).  The  sidereal  mode  is  good  for  astrometry  measurement  and  calibration  but  a 
satellite  may  be  quickly  passing  out  of  the  field  of  view;  this  makes  the  acquisition  of  a  long 
image  sequence  or  long  exposure  time  impossible.  TRM  offers  several  advantages.  First,  the 
satellite  brightness  is  concentrated  within  a  few  pixels,  which  increases  the  signal-to-noise  ratio. 
Second,  the  satellite  remains  almost  at  the  same  image  location,  and  several  images  can  be 
acquired  over  longer  period.  Third,  the  satellite  becomes  the  only  static  object  in  the  image 
sequence,  which  makes  it  easy  to  detect  against  the  moving  background  of  stars.  This  also  makes 
it  easy  to  discriminate  the  flashes  caused  by  cosmic-ray  hits  (CRs).  Flowever,  the  astrometric 
calibration  may  be  less  accurate.  The  star  streaks  must  be  detected  and  their  centroids  calculated. 
The  accuracy  of  the  centroids  calculation  is  limited  by  the  SNR,  which  is  lower  than  the  star  SNR 
in  SSM. 

The  acquisition  of  images  in  TRM  mode  is  more  complex  than  in  SSM  mode.  After  pointing  at 
the  appropriate  location  at  the  appropriate  time,  the  sensor  motion  must  be  adjusted  in  such  a  way 
that  the  residual  motion  rate  produces  the  expected  tracking.  While  tracking,  the  sensor  can 
acquire  several  images,  which  will  produce  more  accurate  and  reliable  measurements  in 
comparison  with  the  acquisition  of  a  single  snap  shot.  Furthermore,  the  TLE  update  process 
requires  a  minimum  of  three  independent  observations;  this  is  a  Space  Surveillance  Network 
‘SSN’  requirement.  Considering  this,  here  is  one  of  the  most  probable  acquisition  scenarios: 


-  slew/settle  to  a  “staging”  point  near  where  the  “tracking  slew”  is  to  begin, 

-  at  the  appropriate  time,  when  the  satellite  enters  the  field  of  view, 

accelerate  &  settle  in  tracking  rate, 

-  3s  exposure  (first  image), 

-  wait  a  few  seconds  before  next  exposure  while  downloading  this  image, 

-  3s  exposure  (second  image), 

-  wait  a  few  seconds  while  downloading, 

-  3s  exposure  (third  image), 

-  wait  a  few  seconds, 

-  3s  exposure  (fourth  image), 

-  slew  to  the  next  satellite  position. 

A  minimum  of  four  images  should  be  acquired.  A  higher  number  would  increase  the  system’s 
sensitivity  and  accuracy.  For  satellites  known  to  be  bright  enough,  three  images  would  be 
adequate  (to  provide  the  three  required  observations)  but  a  fourth  one  is  required  because,  quite 
often,  a  coincidental  event  forces  the  rejection  of  one  of  the  acquired  frame.  This  is  the  case 
when  the  RSO  is  overlapped  by  a  star  streak  or  a  cosmic-ray  hit.  For  fainter  RSOs,  it  would  be 
preferable  to  acquire  a  longer  sequence  (9  to  16  images)  with  a  longer  exposure  times.  This 
increases  the  SNR  (signal  to  noise  ratio  -  by  processing  with  an  image  summation  algorithm)  and 
makes  detection  easier. 

With  a  space-based  system,  jitter  has  a  certain  impact  on  image  quality.  It  is  expected  that  the 
jitter  is  of  the  same  order  as  the  size  of  the  optical  PSF.  In  fact,  the  detection  algorithms  must 
consider  an  effective  PSF,  which  is  the  combination  of  the  optical  PSF  and  jitter.  However,  there 
is  no  intent  to  try  to  separate  the  jitter  effect  from  the  optical  PSF,  and  it  is  not  believed  that  it  is 
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possible.  However,  there  is  a  possibility  that  the  jitter  is  not  constant  over  the  image  sequence, 
and  that  some  images  within  the  sequence  have  a  better  effective  PSF  than  others.  In  such  a  case 
it  would  be  interesting  to  evaluate  the  image  quality  and  reject  the  images  that  are  more  affected 
by  jitter.  However,  this  effect  is  not  considered  in  the  following  algorithm  development.  If  the 
usage  of  a  space-based  platform  proves  that  jitter  variations  are  not  negligible  and  should  be  taken 
into  account,  the  sensor  should  acquire  a  few  extra  images  in  expectation  that  not  all  images  in 
the  sequence  will  be  used  for  detection. 

The  drift  is  also  a  common  source  of  error  that  occurs  during  acquisition.  It  is  caused  by  a 
tracking  rate  that  does  not  exactly  mimic  the  desired  satellite  motion.  It  can  be  taken  into  account 
by  image  processing  algorithms.  First,  it  distorts  the  shape  of  the  PSF.  This  effect  could  be 
corrected  with  ‘sharp’  or  ‘super-resolution’  algorithms,  but  this  is  not  necessary  for  detection 
purposes.  Second,  it  makes  the  satellite  to  appear  at  different  places  in  the  various  images  of  the 
sequence.  This  effect  is  more  important  and  is  taken  into  account  in  the  processing  methodology 
presented  in  the  document.  Its  correction  makes  it  possible  to  apply  a  more  sensitive  global 
detection  algorithm. 
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4  Scene  and  acquisition  models 


It  is  easy  to  conceive  of  ad-hoc  processing  for  the  acquired  image.  However,  this  processing  is 
rarely  optimal.  To  support  a  good  processing  approach,  a  good  comprehension  of  the  image 
formation  process  and  scene  content  is  necessary.  This  chapter  presents  the  image  model  used  for 
the  image  processing  algorithm  development.  However,  before  arriving  at  an  appropriate  image 
model,  it  is  also  necessary  to  model  the  scene  and  other  objects  that  contribute  to  form  the  image 
content. 

The  image  content  has  several  independent  components  that  can  be  processed  separately.  They 
are  the  astronomical  scene,  RSO  element  (which  has  a  different  motion),  cosmic-ray  hits, 
background  illumination  and  sensor  artefacts  (noise  and  bias).  The  complete  image  model  is 
presented  in  Section  5.1.  This  model  may  seem  complex,  but  if  we  look  at  the  individual 
components  first,  then  the  resulting  model  will  seem  more  intuitive. 

4.1  Astronomical  scene  model 

The  astronomical  scene  content  ‘A’  depends  on  the  telescope  field  of  view  ‘FOV’,  camera 
rotation  ‘Rot(t)’  and  the  real  pointing  ‘PA(t)’  at  the  time  tk: 

Ak(Ra,  Dec)  =  Sky  (FOV,  PA(tk),  Rot(tk)) ,  ( 1 ) 

i.e.  Ak  is  a  subset  of  the  global  scene  ‘Sky’.  This  expression  represents  a  star  field  acquired  by 
the  camera.  It  can  also  be  the  result  of  an  image  reconstruction  using  a  star  database  (useful  for 
the  star  recognition  process)  or  an  image  reconstruction  resulting  from  the  processing  of  the  raw 
images. 

If  the  pointing  system  is  moving  and  the  image  is  acquired  over  the  time  period  ‘dt’,  the  resulting 
scene  is  the  cumulative  integration  of  the  signal  over  this  period  of  time,  i.e.  the  interval: 
‘|tk,tk+dt’; 

Ak  (Ra,  Dec  |tk,tk+dt)  =  Sky  (FOV,  PA(tk),  Rot(tk),  |tk,tk+dt )  (2) 

The  static  scene  can  be  seen  as  the  sum  of  all  individual  stars  contained  in  the  FOV: 

A(Ra,  Dec)  =  £  an(Ra,Dec)n  ,  (3) 

where  the  individual  star  ‘n’  has  their  own  intensity  ‘an’  and  position  ‘(Ra,Dec)n’.  For  a  dynamic 
scene,  i.e.,  when  the  telescope  pointing  is  changing  because  of  the  tracking  motion,  jitter  or  drift, 
the  motion  blur  must  be  taken  into  account: 

A(Ra,  Dec)  =  I  an(Ra,Dec)n  0  Z(PA(t)|t,t+dt),  (4) 

where  ‘0’  is  the  convolution  operator  and  ‘Z(PA(t)|t,t+dt)’  is  the  line  segment  (or  other  streak 
shape)  introduced  by  the  pointing  motion  PA(t)  during  the  interval  from  t  to  t+dt. 
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4.2  Image  formation  and  transformation 


Between  the  sidereal  scene  and  the  image  frame,  the  signal  undergoes  several  transformations: 
Global  scene  (RA,Dec), 

->  Selected  scene  (RA,Dec), 

->  Scene  sampling  ( i "  j ”), 

Optical  distortion  (i’,j’), 

->  Other  optical  effect  (i,j):  final  image. 

One  may  note  that  the  optical  distortion  comes  before  the  scene  sampling  and  this  is  true.  But 
these  operations  are  independent  and  their  order  can  be  swapped.  The  advantage  of  this  swapped 
model  is  that  the  geometric  distortion  is  now  a  function  of  the  image  coordinates.  Otherwise,  the 
geometric  distortion  should  be  expresses  in  local  (Ra,Dec)  coordinates  and  this  would  make  the 
representation  more  complicate  because  (Ra,Dec)  are  not  specific  to  the  image.  Pointing  and  FOV 
must  be  taken  into  account.  When  already  in  pixel  coordinates,  the  geometric  distortion  can  be 
easily  corrected  with  a  pixel  interpolation  algorithm  that  uses  a  calibrated  distortion  map  (which 
is  measured  in  laboratory).  In  this  case,  the  scene  sampling  only  means  the  conversion  of  pixel 
index  into  real  space  coordinates.  Similarly,  the  other  optical  effect  (like  the  PSF)  occurs  in  the 
optics  before  the  sampling  by  the  CCD.  But  these  effects  are  also  independent  and  can  be 
modeled  on  the  image  pixel  coordinate  system. 

Accordingly  to  this  model,  when  the  telescope  points  in  a  direction  PA(tk),  it  selects  a  portion  of 
the  sky  ‘Ak’  which  is  delimited  by  the  FOV  and  camera  rotation  ‘Rot(tk)\  This  selects  the  scene: 

Ak(Ra,  Dec)  =  Sky  (FOV,  PA(tk),  Rot(tk))  ( 1 ) 

This  scene  contains  celestial  objects  whose  brightness  is  expressed  in  magnitude  or  radiance 
(W/nr/sr).  For  practical  reasons,  it  is  useful  to  think  in  terms  of  an  equivalent  signal  expressed  in 
pixel  counts  (or  digital  number:  DN).  Hence,  the  intensity  transformation  ‘9K.’  is  a  function  that 
takes  into  account  the  optical  gain  (aperture),  sensor  sensitivity  and  spectral  match  (between  the 
object  and  the  sensor)  and  produces  such  numbers. 


IAk(Ra,  Dec)  =  91  (  Ak  (Ra,  Dec))  (5) 

The  optical  system  creates  an  imperfect  image  of  this  scene,  and  the  image  formation  process  can 
be  conceived  in  three  distinct  steps:  the  transformation  ‘@’  of  the  celestial  scene  into  a  perfect 
image  grid,  the  transformation  ‘O’  of  this  perfect  image  grid  into  a  distorted  image,  and  finally 
the  addition  of  the  optical  aberrations  ‘0’  like  PSF.  Hence,  the  first  step  of  the  image  formation 
is  the  transformation  of  (Ra,Dec)  coordinates  into  a  series  of  temporary  (i",j”)  indices; 


IAk  (i"  j”)  =  @  (IAk(Ra,  Dec), 


(6) 
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where  the  index  ‘k’  indicates  that  this  transformation  is  specific  to  the  image  frame  ‘k’  which  was 
acquired  at  time  ‘tk’.  This  transformation  links  every  image  pixel  to  corresponding  real  sky 
coordinates  (Ra,Dec).  This  image  also  suffers  the  geometric  distortion  ‘O’  caused  by  imperfect 
optics  and  the  resulting  image  is: 

iAk(i’j’)  =  n(iAk(i",j”)).  (7) 

Finally,  the  other  optical  effects  (noise,  PSF,  gain,  etc.)  are  added  within  the  final  transformation: 

IAk(i,j)  =0(IAk(i’,j’)),  (8) 

This  transformation  does  not  affect  the  object  location  like  the  geometric  distortion  ‘O’  and  the 
transformation  ‘@k’  of  the  sidereal  coordinates  into  image  indices.  It  instead  reduces  image 
“sharpness”. 


So,  the  global  transformation  is: 

IAk  (i,j)  =  0  (Q  (  @  (0K  (Ak(Ra,  Dec)) ))), 


(9) 


These  transformations  are  reversible.  Without  considering  the  blur  (caused  by  the  PSF),  the  pixel 
coordinates  are  linked  to  sidereal  coordinates  with: 


IAk(Ra,  Dec)  =  @A  Q1  ( IAk  (i,j)).  (10) 

To  these  transformations,  one  can  add  the  radiometric  transformation  ‘011’  that  converts  the 
object’s  brightness  (initially  expressed  in  magnitude)  into  pixel  count,  but  this  transformation  is 
not  required  in  this  document  because  all  image  processing  algorithms  only  use  the  pixel  count 
notation.  Let  us  consider  that  this  brightness  to  pixel-count  transformation  function  is  integrated 
somewhere  within  the  optical  function  ‘0’. 

4.3  Astronomical  calibration  of  an  image 

There  are  several  commercial  software  packages  like  ‘CCDSoft’  or  ‘PinPoint’  (Refs,  5  and  6)  that 
can  detect  stars  within  an  image,  compare  with  a  star  database  and  evaluate  image  features.  For 
example: 

Software:  ‘PinPoint’ 

Input:  corrected  image  =  f>  '(IA) 

Outputs:  FOV,  PA(t),  Rot(t)  and  transformation  method  @k. 

For  practical  reasons,  let  us  call  this  operation  the  function  ‘Astrometry’: 

[  @k,  FOV,  PA(t),  Rot(t)  ]  <r  Astrometry  (Q  _1(  IAk  ) ).  (1 1) 

This  notation  will  be  used  later.  For  an  accurate  evaluation  of  the  astrometric  measurements,  note 
that  the  image  must  be  previously  rectified  (to  remove  geometric  distortions). 
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4.4  Object  position  model 


Accordingly  to  the  TLE  and  prediction  model,  the  predicted  apparent  position  of  the  observed 
RSO  (designated  by  ‘R’)  is  the  vector: 

Rp(TLE,t)  =  (Ra(t),  Dec(t)  )RP  ,  (1 2) 


where  ‘Rp’  is  the  vector  and  ‘(Ra,  Dec,  t)Rp’  the  elements  of  this  vector.  However,  because  of 
orbit  perturbations,  the  real  RSO  position  is: 

RA(t)  =  (Ra(t),  Dec(t)  )RA  (13) 

The  uppercase  index  ‘p’  means  ‘predicted’  and  ‘A’  means  an  observation  astrometrically 
calibrated.  The  ‘Ra’  and  ‘Dec’  coordinates  of  the  RSO  are  expressed  as  a  function  of  time 
because  the  object  and  sensor  are  moving. 

Figure  2  shows  an  example  where  the  real  position  is  not  exactly  the  predicted  one.  It  also 
indicates  that  the  sensor  is  not  able  to  point  exactly  at  the  desired  position,  because  of  limited 
maneuvering  capabilities.  In  addition,  the  sensor  does  not  point  exactly  at  what  is  believed, 
because  of  limited  sensor  accuracy. 


Figure  2:  Estimated  and  real  object  positions  and  pointing  vectors 


4.5  Object  image  model 

Similarly  to  the  astronomical  content,  the  RSO  can  be  projected  to  the  image  plane  with 
transformation  functions  ‘O’  and  ‘@’  and  becomes  the  image  component  ‘IR\  Its  absolute 
position  is  ‘RA(tk)’  while  its  apparent  position  in  the  image  is  relative  to  the  telescope  pointing, 
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i.e.;  FR(t)  =  RA(t)  -  PA(t).  It  also  has  an  intensity  ‘aR’,  which  generates  an  equivalent  signal  upon 
the  image  frame.  The  astrometric  conversion  function  @k  converts  the  RSO  sidereal  coordinates 
directly  into  pixel  indices.  If  the  relative  position  ‘FR’  is  used,  then  the  pointing  value  must  be 
included  in  the  conversion. 


IR(tk)  =  0tc/aR)  0  £2  @k  (  RA(tk) )  ( 1 4A) 

or: 

IR(tk)  =  01t(aR)  0  Q  @k  [(RA(tk)) )  0  Z(tk,tk+dt)  ]  (14B) 

‘aR’  is  the  RSO  magnitude  while  ‘9lt(aR)’  is  its  equivalent  in  pixel  counts.  The  function 
‘Z(tk,tk+dt)’  represents  the  shape  of  the  RSO  streak  caused  by  its  apparent  motion.  In  tracking 
mode,  it  is  assumed  to  be  Dirac’s  delta  because  ‘PA’  is  supposed  to  be  equal  to  ‘RA’.  In  sidereal 
mode  this  shape  depends  only  on  ‘d/dt  RA(tk)’  and  should  be  a  line  segment.  However,  the  jitter 
and  drift  may  make  this  function  more  complex.  Its  exact  representation  is: 

Z(tk,tk+dt)  =  d/dt  [RA(tk)-PA(t)]  dt  =  t  d/dt  FR(t)  dt  (15) 

4.6  Pointing  model 


Theoretically,  the  telescope  should  be  maneuvered  to  point  at  the  RSO  at  its  predicted  position, 
i.e.  the  expected  pointing  is: 

PR(tk)  =  Rp(TLE,tk)  =  (Ra(tk),  Dec(tk)  )RP.  ( 1 6) 


In  practice,  the  pointing  system  has  some  inaccuracies  but  can  achieve  a  adequate  pointing  to 
ensure  that  the  expected  RSO  is  close  to  the  optical  axis.  Using  the  attitude  measurement  systems 
(magnetometers,  accelerometers,  star  tracker,  etc.),  the  real  pointing  ‘P0,  parameters  are  measured 
and  recorded: 


P°(tk)  =  (Ra(tk),  Dec(tk))Po  ~  (Ra(tk),  Dec(tk)  )RP 


(17) 


This  pointing  is  only  an  estimation  and  can  have  a  relatively  large  error  component.  The  upper 
script  ‘°’  means  ‘Original’  pointing  estimation  while  ‘R’  means  desired  pointing  (at  the  RSO). 
However,  after  an  image  is  acquired,  the  image  content  is  checked  and  the  known  stars  are 
identified  (using  comparisons  with  a  database  of  reference).  The  real  pointing  parameter  can  be 
evaluated  with  almost  the  same  accuracy  as  the  star  position  accuracy  in  the  catalogue. 

PA(tk)  =  (Ra(tk),  Dec(tk)  )PA  (18) 


This  pointing  estimation  (with  the  upper  script  ‘A’)  is  astrometrically  calibrated  and  it  is  deduced 
from  the  image  content  with  software  like  ‘Pinpoint’,  as  described  previously  in  Eq.  11. 


DRDC  Valcartier  TR  2009-095 


11 


4.7  Tracking  modes 


In  sidereal  tracking  mode,  the  telescope  points  at  the  object  and  attempts  to  stop  its  apparent 
sidereal  motion. 

P°(tk)  =  (Ra,  DecjpO.  (19) 

In  RSO  tracking  mode,  the  telescope  keeps  moving,  trying  to  maintain  the  predicted  RSO 
position  ‘Rp’  in  the  middle  of  the  FOV.  For  a  short  period  of  time,  the  telescope  slew  motion  is 
set  to  constant  ‘Ra’  and  ‘Dec’  rates  (TRM:  track  rate  mode). 

P°(tk)  =  (Ra(tk),  Dec(tk))P0  «  Rp(TLE,t)  =  (Ra(tk),  Dec(tk))Rp.  (20) 

For  long  acquisition,  active  motion  sensing  and  tracking  error  correction  should  be  applied,  but 
this  is  not  considered  in  the  current  case  where  the  total  acquisition  time  is  not  long  enough.  The 
resulting  drift  is  not  negligible  but  it  is  taken  into  account  and  compensated  in  the  processing 
described  below. 

4.8  Pointing  offset  and  apparent  position 

The  pointing  offset  ‘Fp(t)’  is  error  which  occurs  between  the  real  pointing  ‘PA(t)’  (determined 
based  upon  astrometry)  and  the  desired  pointing  which  is  the  predicted  RSO  position  ‘Rp(TLE,t)’. 
It  is  obtained  with  the  vector  subtraction; 


Fp(t)  =  Rp(TLE,t)  -  PA(t),  (2 1 ) 

where  ‘PA’  and  ‘RP’  are  unit  vector  but  not  ‘Fp’. 

Similarly,  the  apparent  position  of  the  RSO  is  the  difference  between  its  real  position  and  the  real 
telescope  pointing  parameters: 

FR(t)  =  RA(t)  -  PA(t)  (22) 

This  is  the  value  the  data  acquisition  and  processing  will  determine. 

4.9  Drift  estimation 

Drift  ‘D’  appears  when  the  pointing  is  changing  during  the  acquisition  because  of  residual 
unpredicted  pointing  motion.  An  example  is  illustrated  in  Fig.  3.  In  sidereal  mode,  the  camera 
should  stare  at  the  stars.  The  pointing  should  be  constant  from  one  frame  to  another.  Usually, 
star  trackers  are  used  and  the  pointing  accuracy  is  maintained  with  a  dynamic  feedback  loop. 
However,  in  TRM,  the  motion  is  estimated  and  the  satellite  is  maneuvering  consequently.  It  is 
virtually  impossible  to  verify  the  tracking,  or  at  least  not  with  an  accuracy  similar  to  that  provided 
by  star  trackers.  With  the  astrometric  calibration,  the  real  pointing  position  can  be  measured  and 
the  pointing  motion  can  be  evaluated  by  comparing  the  pointing  of  different  frames: 


D(t)  =  dPA(t)/dt  dt  or  dPA(t)/dt  t 


(23A) 
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or 


D  (tk)  =  PA(tk)  -  PA(tk_0  . 


(23B) 


where  ‘k’  is  the  image  frame  index  and  ‘tk’  the  time  of  acquisition  (not  the  exposure  time)  and  dt 
the  integration  period  (which  is  the  exposure  time  Y).  This  drift  value  ‘D(tk)’  will  be  used  later 
in  the  image  processing  section  for  the  realignment  of  images. 

In  TRM,  a  similar  drift  can  be  estimated.  It  is  assumed  that  the  RSO  position  is  only  known 
approximately,  but  its  motion  (relative  to  the  sensor)  is  known  with  accuracy.  Thus,  assuming 
that  the  camera  is  not  rotating,  the  drift  can  be  estimated  by  comparing  the  real  pointing  relative 
to  the  estimated  RSO  position  (i.e.  the  pointing  offset)  in  time; 

D(t)  =  d/dt  -FR(t)  dt  =  d/dt-  Fp(t)  dt  =  d  (-  PA(t)  +  Rp(TLE,t))  /dt  dt  (24A) 

or: 

D(tk)  =  [  -PA(tk)  +  PA(tk_0  +  Rp(TLE,tk)  -  RP(TLE,tk.|)  ]  (24B) 


Once  the  astrometric  conversion  ‘@k’  is  known,  the  drift  can  be  translated  into  pixel  (i,j) 
displacement  per  frame,  i.e.:  (di/dk,dj/dk).  Because  the  optical  axis  is  the  image  central  pixel 
(io,jo),  this  can  also  be  thought  of  as  corresponding  to  the  pointing  vector  translated  into  image 
coordinates: 


(io,jo)  -  O  @k  ( PA(tk), 

the  relative  displacement  provided  by  the  drift  is: 

Dk  =  (di/dk,dj/dk)  =  [  Q  @k(PA(tk)  +  D(tk) )  ]  -  (i0,j0). 


(25) 


(26) 


In  summary,  it  is  important  to  note  that  this  drift  estimation  method  has  the  ability  to  compensate 
for  two  important  sources  of  error.  First,  for  cases  where  the  optical  system  can  only  achieve 
poor  pointing  accuracy,  this  does  not  affect  the  processing  capability,  since  the  drift  is  estimated 
using  the  astronomical  image  content.  Second,  if  the  RSO  motion  has  important  changes 
(because  of  the  relative  motion  between  the  RSO  and  the  observing  telescope),  this  is 
automatically  taken  into  account  by  propagation  model  (which  calculates  RP(TLE,t)).  Hence, 
this  set  of  equations  facilitates  evaluating  the  exact  drift  and  relative  displacement  and  makes 
possible  exact  image  realignment. 

4.9.1  Detecting  several  co-located  RSOs 

Usually,  co-located  RSOs  have  almost  the  same  orbital  parameters.  For  a  short  period  of  time, 
they  do  not  have  apparent  relative  motion.  They  move  as  a  single  body.  Hence,  the  previous  drift 
estimation  will  allow  the  detection  of  such  a  group  of  RSOs  (common  in  GEO  belt).  However,  if 
the  apparent  relative  motion  becomes  apparent,  the  previous  set  of  equation  can  be  adapted 
simply  by  defining  a  position  vector  ‘Rp(TLE,t)’  for  each  one  the  expected  RSO.  This  will  be 
more  CPU  intensive  but  the  case  remains  manageable. 
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Drift-compensated  realigned  images 


Figure  3:  Drift  measurement,  in  tracking  mode,  relatively  to  the  expected  RSO  positions 
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5  Image  model  and  processing 


The  previous  chapter  provides  us  with  information  that  can  be  used  in  image  processing 
algorithms;  the  scene  and  RSO  image  components.  It  remains  to  add  the  none-image  components 
to  the  image  model:  sensor  noise  and  bias,  cosmic-ray  hits  and  diffuse  light  sources.  With  a 
complete  model,  it  will  be  possible  to  establish  a  method  for  the  extraction  of  a  specific 
component  such  as  the  RSO,  which  is  the  goal  of  this  study. 

5.1  The  image  model 

The  image  content  can  be  described  using  the  following  model.  The  original  image  acquired  by 
the  sensor  (with  the  superscript  ‘°’  for  ‘original’)  depends  on  the  telescope  pointing  parameters 
‘PA(tk)’  at  the  time  ‘tk’,  the  CCD  temperature  ‘T’,  the  exposure  time  7’  (which  replaces  the  time 
interval  dt),  the  FOV  and  other  optical  parameters  hidden  inside  the  optical  function  ‘0’.  For  an 
image  sequence,  the  ‘kth’  images  of  the  sequence,  acquired  at  time  ‘tk’,  is: 

I°k  =  I°(T,tk7,  PA(tk)))  =  IA(tk,t)  +  IR(tk,t)  +  Bk(T,t,  PA(tk))  +  nk(T,t)  +  Ck(t). 

(27) 

‘I0,  is  the  summation  of  the  independent  astronomical  ‘IA  ‘  and  RSO  ‘IR  scenes,  background  ‘B’, 
noise  ‘n’  and  cosmic -ray  hits  ‘C’.  The  astronomical  and  RSO  content  7A  and  ‘IR’  images  are 
previously  defined  in  Eqs.  9  and  14.  This  image  model  should  also  include  a  description  of  hot 
and  dead  pixels.  In  practice,  this  information  can  be  recorded  on  static  maps  and  the  images  can 
be  corrected  in  a  pre-processing  step. 

The  goal  of  the  following  image  processing  techniques  are  to  separate  these  independent 
components  and  isolate  the  IR(tk,t)  component,  which  is  necessary  for  achieving  detection. 

5.2  Background  estimation  and  removal 

The  image  background  ‘B’  is  principally  caused  by  sensor  artifacts  ‘Dk(T,t)\  and  depend  on 
temperature,  exposure  time  and  other  sensor  characteristics  (like  the  dark  current,  DC  offset, 
ADC  noise,  etc.).  However,  it  may  also  contain  external  sources  ‘Lf  like  zodiacal  light  (or 
diffuse  nebula)  or  stray  light  ‘L2’  and  these  components  change  with  the  telescope  pointing. 

Bk  (T ,t,  PA(t))  =  2>k(T,0  +  L1(PA(tk),  f)  +  L2(PA(tk),  t)  (28) 

Stray  light  can  be  produced  by  strong  light  sources  close  to  the  field  of  view.  Particular  effort  is 
applied  to  stay  away  from  the  Sun,  the  Moon  and  the  Earth’s  horizon.  Very  bright  stars  like  Vega 
or  Sirius  can  also  produce  visible  stray  light  and  significantly  contribute  to  the  background.  With 
a  very  good  optical  design  and  meticulous  operation,  L|  and  L2  should  not  be  significant  light 
components  and  can  be  ignored  most  of  the  time,  but  there  is  no  guarantee  they  will  never  appear. 


This  background  ‘Bk’  is  harmful  to  the  image  compression  or  other  image  processing  algorithms. 
‘Bk’  must  be  estimated  (‘<Bk>’)  and  removed.  Good  background  estimation  methods  are 
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presented  in  Ref.  4.  A  ‘good’  estimation  method  is  one  that  produces  a  minimal  residue  ‘3t’,  i.e., 
the  difference  between  the  real  background  and  the  estimation,  is  less  than  the  noise  amplitude; 


Stk  =  |  Bk  (T,t,  PA(tk))  -  <Bk  (T ,t,  PA(tk))>  |  <  cn,  (29) 

There  are  several  methods  to  evaluate  image  background  ‘<B>’  (Refs.  1  to  4).  The  most  popular 
remains  the  measurement  of  the  dark  frame  ‘Dk(T,t)’  (Ref.  3)  by  acquiring  an  image  with  a  closed 
shutter.  However,  if  Li  or  L2  is  present,  it  is  not  included  in  this  background  estimation.  Even  if 
Li  and  L2  are  null,  there  is  no  guarantee  that  the  dark  frame  subtraction  method  will  leave  a 
residue  at  an  acceptable  level.  Reference  3  demonstrated  that  dark  frames  acquired  with  a  very 
small  temperature  difference  may  leave  a  large  residual  dark-frame  count,  i.e.; 


Slk  =  |  5)(Tk,t)  -  5)(Tk+dT,t)  |  >  on  if  |  dT  |  >  ‘temperature  tolerance’  (30) 

With  the  tested  CCD  (in  Ref.  3),  the  temperature  tolerance  (to  obtain  an  acceptable  residue)  was 
no  more  that  0.1C  and  this  was  less  than  the  temperature  measurement  accuracy.  This  was  an 
unsatisfactory  CCD.  Note  that  there  is  always  a  small  temperature  fluctuation  caused  by  the 
feedback-loop  hysteresis  between  the  temperature  measurement  and  cooler-system  control.  It  is 
possible  that  the  measured  temperature  be  never  exactly  the  CCD  temperature  (depending  on  the 
location  of  the  temperature  sensor). 

A  possible  solution  to  minimize  the  residue  is  the  acquisition  of  several  dark  frames  (at  different 
temperatures  ‘Tn’)  combined  with  a  very  accurate  CCD  temperature  measurement  system.  The 
appropriate  dark  frame  could  be  interpolated  from  these  references; 


5>k(Tk,0=  a  S>  (T,  ,t)  +  b  ®(T2,0, 

(31  A) 

b  =  (Tt-T1)/(T2-T1),  and 

(3  IB) 

a  =  1-b. 

(3 1C) 

This  method  is  the  most  economical  in  terms  of  CPU  capability  and  may  be  used  onboard  a 
satellite.  Usually,  the  same  dark  frame  is  used  for  an  entire  image  sequence,  but  if  temperature 
changes  during  the  sequence  acquisition,  there  is  a  possibility  of  adapting  the  dark  frames  for 
every  individual  image  frame. 

The  other  methods  (Refs.  1  and  4)  consist  of  evaluating  the  background  by  software  which 
assumes  that  the  background  is  a  smooth  function  while  the  celestial  objects  are  sharp  and  always 
brighter  than  the  background.  If  Li  and  L2  are  present  and  correspond  to  this  background 
characteristic,  they  are  automatically  included  in  the  background  estimation.  References.  1  and  4 
describe  iterative  methods  based  on  polynomial  fit  or  local  statistics.  These  methods  converge 
quickly  toward  very  accurate  estimations.  The  background  residue  produced  by  these  methods  is 
less  than  0.2  an.  However,  they  are  CPU  intensive,  and  sufficient  computing  capacity  may  not  be 
available  on  a  satellite.  However,  they  are  adequate  for  deferred  processing  performed  by  ground 
stations. 

The  current  case  (image  sequence)  has  another  characteristic  that  is  neither  exploited  by  the  dark 
frame  method  nor  the  iterative  local  estimation  of  Refs  1  and  4.  This  static  background  is 
systematically  present  in  every  frame  of  the  image  sequence.  In  TRM  mode,  the  background  is 
almost  the  same  in  every  frame,  and  the  stars  never  appear  at  the  same  place.  Only  the  satellite 
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appears  at  the  same  place  if  it  is  perfectly  tracked.  A  first  background  estimation  can  be  obtained 
by  extracting  the  minimum  of  the  image  sequence: 

Imin  =  minimum  (I°i,  . . . ,  I°K  )  (32) 

The  drawback  of  this  estimation  is  that  it  underestimates  the  background  because  of  a  negative 
bias  introduced  by  the  systematic  selection  of  all  negative  noise  peaks.  It  also  includes  all  static 
objects  like  very  well  tracked  RSO.  A  median  filter  should  produce  better  results; 

Imed  =  smoothing  (median  (I°i,  . . .,  I°K  ))  (33) 

This  estimation  should  ignore  single  bright  events  like  the  presence  of  star  streaks  or  cosmic-ray 
hits  in  a  single  frame.  A  smoothing  operation  may  be  required  to  obtain  a  background  without 
texture.  Equation  31  should  also  provide  a  better  background  average  because  it  does  not 
systematically  extract  the  minimum  noise  spikes.  However,  it  will  always  contain  static  objects. 
So,  the  smooth  background  assumption  must  be  used  again,  i.e.  a  static  object,  like  a  well-tracked 
RSO,  is  small  (compared  with  the  image  size),  relatively  sharp  (almost  the  PSF  shape)  and 
brighter  than  the  background  while  the  background  is  smooth  and  dark.  So,  the  following 
assumptions  can  be  made: 

1)  a  smoothed  background  still  represents  the  background, 

2)  a  smoothed  image  (RSOs  +  Background)  contains  a  background  with  a  very 
attenuated  RSOs, 

3)  the  minimum  between  the  background  and  the  smoothed  background  remains  the 
background,  except  for  point  objects  (RSO)  which  are  severely  attenuated. 

The  background  estimate  is  obtained  with  the  function; 

<Bk>  =  minimum  (Imed  ,  Imed  0  Targe  impulse  response’).  (34) 

By  repeating  the  process  with  an  iterative  scheme,  the  RSO  completely  disappears,  leaving  only 
the  background.  This  method  is  described  in  detail  in  Ref.  4  where  the  background  is  evaluated 
in  the  presence  of  stars  and  where  the  stars  are  eliminated  from  the  background  by  using  this 
approach: 


<Bk[>  =  Imed,  (i.e.  the  first  estimate)  (35) 

<Bkn>  =  minimum  (<Bkn_i>  ,  <Bkn_i>  0  Targe  impulse  response’).  (36) 

Once  the  background  is  adequately  estimated,  a  background-free  image  TB’  can  be  produced 
(Ref.  1  to  4). 

l\  =  I°k  -  <Bk>  (37) 

TBk  =  IAk  +  IRk  +  nk  +C k  +  5£k  (38) 

and  ~  0.  The  upper  script  ‘_B’  means  background  removed. 
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5.3  Star  detection 


The  next  step  consists  of  detecting  the  stars.  This  detection  will  have  several  usages.  First,  the 
detected  stars  are  used  to  measure  astrometry  parameters,  i.e.,  the  real  pointing  values  of  the 
telescope  PA(tk)  along  with  the  measurement  of  the  image  rotation.  Second,  astrometry  is  used  to 
measure  the  drift  (in  satellite  tracking  mode)  of  the  telescope  pointing  relative  to  the  expected 
satellite  position.  Finally,  the  presence  of  stars  is  a  “determinate  false  alarm  source”  relative  to 
satellite  detection.  The  better  we  are  able  to  predict  and  subsequently  remove  false  alarms,  the 
easier  RSO  detection  will  be. 

Depending  on  the  tracking  mode,  the  stars  appear  like  streaks  or  points  (with  a  PSF  profde).  In 
Ref.  1,  methods  were  developed  for  the  detection  of  these  two  kinds  of  objects.  Sidereal  tracking 
mode  is  well  solved  in  Ref.  1.  Here,  the  effort  will  be  focused  on  constant  track  rate  mode 
(TRM). 

In  tracking  mode,  the  stars  appear  in  the  image  as  streaks.  They  also  appear  fainter  because  their 
energy  is  spread  over  more  pixels.  Thus,  the  detection  methods  using  clipped  signal  and  Hough 
transform  may  not  be  appropriate  when  there  are  not  enough  bright  stars  in  the  FOV.  In  such  a 
case,  the  matched  fdter  described  in  Ref.  1  is  more  appropriate  for  this  task.  This  method  is  able 
to  detect  a  single  streak  as  faint  as  ISNR=0.5,  i.e.,  an  average  signal-to-noise  ratio  of  0.5  per 
pixel.  However,  such  a  sensitivity  may  not  be  required  if  there  are  enough  very  bright  stars. 

The  detection  method  of  Ref.  1  can  produce:  1)  an  image  T+s’  that  contains  only  the  star  streaks, 
2)  the  list  of  the  centroids  for  each  one  of  these  streaks  ‘Sk(n,i,j)’,  and  3)  a  binary  image  mask 
‘Msk’  which  locates  the  streak  pixels  in  the  background-free  image  TBk\ 

Streak_detector(FBk)  ->  [I+Sk,  Sk(n,i,j),  Msk  ]  (39) 

The  follwing  describes  the  method  explained  in  Ref.  1.  The  background-free  image  TBk’  is 
convolved  with  the  expected  streak  shape  ‘Z(PA(tk)|tk,+/)\  The  convolved  image  ‘Ak’  contains 
convolution  peaks  which  indicate  where  the  streaks  are. 

Ak  =  FBk  ®  Z(PA(t)|t,t+dt)  (40) 

The  iterative  convolution  and  clipping  method  of  Ref.  1  allows  improving  filtering  performance 
(‘m’  is  the  iteration  index)  ; 

Akm  =  FBkm  0  Z(PA(t)|t,t+dt),  (41) 

and 

FBk,m+i)  =  minimum  {I‘Bkm,  2  Akm  }.  (42) 

Once  the  iterative  process  converges  toward  a  stable  solution  (about  three  iterations  normally 
improves  significantly  the  performance),  the  final  convolved  image  is  clipped  at  a  predetermined 
threshold  to  produce  a  binary  extraction  mask  ‘Msk\  In  the  case  of  multiple  streak  detection,  this 
threshold  should  be  a  value  close  to  the  noise  level  ‘an’; 

Msk  =  Ak  >  an.  (43) 

The  image  T+Sk’  containing  only  the  star  streak  is  obtained  by  multiplying  this  mask  by  the 
background-free  image: 
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I+Sk  =  MSk  *  rB.  k  (44) 

This  image  should  contain  all  major  astronomical  objects  and  is  a  very  good  estimation  of  the 
astronomical  scene  defined  in  the  image  model,  i.e. : 

IA(tk,0  ~  I+Sk  .  (45) 

For  the  localization  of  every  individual  star,  the  mask  ‘Msk’  is  segmented  into  ‘N’  individual 
objects  ‘Mskn’.  The  segmentation  method  is  presented  in  Ref.  2; 

Msk  =  SMskn.  (46) 

For  each  streak,  a  sub-image  T+Skn’  is  extracted  by  multiplying  the  background-free  image  by  the 
corresponding  single-object  mask; 

I+Skn  =  Mskn  *  Fb  (47) 

The  center  of  intensity  ‘i,j’  (same  calculation  as  the  center  of  mass)  and  total  intensity  ‘an’  is 
calculated  for  every  ‘n’  streaks  and  a  global  list  ‘Sk(n,i,j,an)’  is  produced; 


Sk(n,i,j,an)  =  [n,  (i,j)=center  of  mass  { I+Skn  },  an  =  total  intensity  { I+Skn  }  ].  (48) 


5.4  Star  streak  removal 


Once  the  star  streaks  are  detected  in  the  image  and  the  mask  ‘Msk’  produced,  the  image  without 
stars  is  obtained  via: 

rsk=rBk(i-  m\).  (49) 

The  super  script  ‘‘s’  means  the  streaks  are  removed  from  the  image  (the  annotation  TB‘sk’  could 
also  be  used).  According  to  the  image  model,  this  image  should  contain  only  the  RSO  plus  the 
noise  and  cosmic-ray  hits; 

I  Sk  ~  IR(tk,0  +  nk(T,Z)  +  Ck(t).  (50) 

Flowever,  Eq.  47  may  leave  signal  artifacts.  When  the  streaks  have  extended  PSF  profiles  and  the 
mask  is  too  narrow,  only  the  middle  part  of  the  streak  is  erased,  leaving  truncated  signals  with 
halo  shapes  around  the  erased  objects.  With  the  multi-frame  detection  method  explained  in 
Section  5.6,  these  artifacts  should  not  have  serious  impacts  on  the  detection  results.  But  if  this  is 
the  case,  Ref.  1  indicates  a  more  accurate  method  for  the  removal  of  PSF-shaped  objects  that  do 
not  leave  such  artifacts.  Removing  such  PSF  shaped  objects  will  require  additional  calculations, 
but  it  will  not  present  a  technical  challenge. 

The  removal  of  the  star  streak  should  reduce  the  production  of  false  alarm,  later  in  the  RSO 
detection  phase.  Flowever,  it  is  not  sure  that  this  is  necessary;  the  RSO  detector  can  be  design  to 
detect  the  RSO,  even  in  presence  of  stars,  as  explained  later  in  section  5.7.  This  is  a  question  of 
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performance,  evaluated  by  testing  with  simulated  data  first  and  then  with  real  data,  that  will 
determine  is  the  removal  of  the  star  streaks  increase  the  detection  reliability. 

5.5  Image  alignment  with  drift  correction 

After  the  removal  of  star  streaks,  it  should  be  relatively  easy  to  detect  the  RSO  in  the  residual 
image.  However,  the  cosmic  rays  produce  strong  artefacts,  usually  brighter  that  faint  RSO.  The 
method  for  discriminating  between  the  RSO  and  cosmic  rays  consists  of  analyzing  the  entire 
image  sequence.  In  tracking  mode,  the  RSOs  should  systematically  appear  at  the  same  place 
while  the  cosmic-ray  hits  appear  at  random  places.  The  RSO  steady-state  condition  is  respected 
only  in  the  absence  of  drift  or  with  compensated  drift.  Thus,  in  the  presence  of  a  drift,  the  image 
sequence  needs  to  be  realigned  before  comparing  the  pixels. 

With  the  previously  defined  drift  (the  difference  between  the  real  pointing  PA(t)  and  the  RSO 
expected  position  Rp(TLE,t) ),  i.e. : 

Dk  =  (di/dk,dj/dk)  =  [  Q  @k  (PA(tk)  +  D(tk) )  ]  -  (io jo).  (26) 

it  is  possible  to  calculate  the  cumulative  drift  ‘X’  (from  first  to  last  frame)  and  realign  the  image 
sequence  frames  relative  to  the  first  frame  (the  ‘X’  mnemonic  is  used  due  to  running  out  of  letters 
in  the  alphabet,  sorry  for  this  bad  mnemonic): 

k 

xk  =  (Xik,  Xjk)=  ^ Dk ,  where  1  <  k  <  K,  (51) 

i 

where  ‘K’  is  the  total  number  of  images  in  the  sequence. 

The  image  frame  Tsk’  can  be  realigned  by  shifting  the  image  by  an  equivalent  cumulative  drift 
displacement; 

I'Dk  =  I'sk  (i-  Xik,  j-Xjk).  (52) 

Hence,  all  the  realigned  images  of  the  sequence  TDk’  should  present  the  RSO  at  the  same  pixel 
location.  To  make  a  long  story  short,  the  upper  case  ‘‘D’  means  the  drift  is  now  removed  or 
compensated.  For  the  complete  story,  one  should  write:  TB's'Dk’  because  the  background  and 
stars  were  also  previously  removed,  but  this  would  only  produce  a  more  confusing  notation. 

The  interest  of  the  image  alignment  is  that  all  the  image  frames  can  be  summed  to  produce  a 
single  image  Ts’  where  the  RSO’s  SNR  is  improved  by  a  factor  ‘a/K’.  This  assumes  a  Gaussian 
noise,  which  is  very  close  to  reality. 

K 

I2  =  Z  I_Dk  (53) 

k= 1 


According  to  the  image  model  of  Eq.  27,  which  is  reduced  to  Eq.  50  after  the  removal  of 
background  and  stars,  the  sum  of  the  realigned  images  is: 
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(54) 


K 

f^Y  [  IR(tk,0  +  nk(T,0  +  Ck(0  ]  . 

k= 1 

K 

~  K  IR(tk,0  +  ns  (T ,t)  +  Y  Ck(0-  (55) 

k=\ 

Using  the  noise  addition  rule  where  a22  =  V(  a  i 2  +  a  2)',  G„y  ~  Vk  onk.  The  signal  amplitude 
increases  with  the  number  of  sample  images  ‘k’,  while  the  noise  amplitude  only  increases  with 
the  Vk  rule.  For  example,  with  an  original  RSO  SNR  =  4,  which  is  very  border-line  for  the 
detection,  the  new  SNR  is  improved  to  ‘12’  for  a  sequence  of  nine  images,  which  results  in  a 
more  readily  detected  RSO. 

5.6  Detection  of  the  RSO 

With  this  last  image,  the  challenge  is  the  detection  of  a  signal  above  the  noise  level  with  the 
rejection  of  the  often  brighter  cosmic-ray  hits.  Let  us  consider  the  following  signal 
characteristics: 


The  RSO  is  almost  constant  in  every  image.  Only  the  possible  detections  are 
considered,  i.e.,  RSO  brighter  than  the  noise  level. 

Cosmic-ray  hits  can  be  very  bright  but  appear  in  the  image  sequence  at  random. 
Lienee,  the  first  detection  is  perfomed  using  the  summed  image  with  a  high  detection  threshold 


Aij  =  I2(i,j)  >  hh  .  (56) 

‘A’  is  a  binary  alarm  matrix.  Typically,  hh  is  greater  than  2  to  3  times  the  noise  level,  which  is 
VKct„  in  the  summed  image.  So,  a  typical  hh  is: 

hh  =  3VKg„  (57) 

After  the  detection  with  the  summed  image,  the  detected  object  can  be  track  back  in  individual 
image  frames  and  detected  again  with  a  lower  detection  threshold  value  ‘hf,  which  is  applied 
only  within  the  region  of  interest  identified  by  the  global  detection.  Lienee,  for  a  global  detection 
‘Ay’,  there  is  a  series  of  individual  local  detections  ‘Aljk\  Hence,  the  complete  detection  test  is: 

Aij  =  [  IE(i,j)  >  hh  ]  and  [  n  (rDk(i,j)  >  h,)  ],  (58) 

k= 1 


where  ‘IT  is  the  product  of  all  logical  conditions,  i.e.  equivalent  to  a  ‘AND’  logical  operator. 

For  a  sequence  of  nine  images,  the  ‘hh’  could  be  9an  while  ‘hf  could  be  as  low  as  2an  (even  lan 
should  be  tested).  The  exact  threshold  value  will  have  to  be  fine  tuned  with  real  data  and  false 
alarm  tolerance  constraints.  However,  9an  and  2an  are  good  assumptions  and  they  will  be 
considered  as  values  of  reference. 

The  limit  of  sensitivity  of  this  filter  is  SNR>3  (when  9  images  are  summed,  the  integrated  SNR 
triples  and  pass  the  test  of  Eq.  58).  It  can  reject  cosmic-ray  hits  of  any  intensity.  Table  I  indicates 
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the  behaviour  of  the  combined  tests.  To  detect  a  RSO,  the  global  and  all  individual  tests  must 
succeed.  The  very  faint  (SNR>3)  RSO  will  pass  the  global  test  and  all  tests  in  the  individual 
frame.  Once  in  a  while,  the  noise  will  succeed  the  detection  test  in  a  frame  (and  generate  a  false 
alarm)  but  never  the  global  test  in  the  summed  image.  A  cosmic-ray  hit  may  be  so  intense  that  it 
may  succeed  the  global  test,  but  in  the  second  test,  most  of  the  individual  frames  (which  contain 
only  noise)  will  fail.  Perhaps,  once  in  a  while,  a  combination  of  noise  and  background  residue 
will  pass  the  tests.  In  such  a  case,  the  background  residue  probably  cover  an  extended  area  and 
the  non-PSF  shape  can  be  used  to  reject  similar  cases. 

One  may  be  tempted  to  use  only  the  second  test,  but  the  combination  of  both  tests  makes 
detection  more  reliable  and  more  sensitive.  The  noise  and  background  residue  can  create  false 
alarms  with  the  second  test  but  will  never  trigger  the  first  one.  With  a  higher  threshold  value  ‘hi’ 
(which  would  reduce  the  false  alarm  rate),  the  combination  of  the  faint  RSO  signal  and  noise  may 
fail  the  second  test.  This  would  make  the  system  less  sensitive  and  this  is  not  desired.  Therefore, 
the  combination  of  both  tests  assures  sensitivity  and  false  alarm  rejection. 

5.7  Detection  of  a  RSO  in  presence  of  a  star  streak 

It  is  possible  that  the  RSO  overlaps  a  star  streak  in  one  (or  more)  of  the  frames.  If  this  occurs, 
this  RSO  signal  should  be  erased  with  the  streak  with  the  consequence  that  the  detection  with  Eq. 
58  fails.  This  situation  can  be  corrected  by  inclusion  of  the  star  streak  mask  within  the  detection 
algorithm.  First,  the  streak  mask  ‘Msk’  must  be  aligned  with  the  realigned  image: 

Ms'Dk  =  Msk(i-Xlk,j-Xjk).  (59) 

With  this  modified  detection  rule; 


Aij=  [  I  (i,j)  >hh] 


K 

and  [  ]^[ 


{  (I  Dk(i,j)  >  hi)  }  or  {  Ms"°k(i,j)  =1  }  ], 


k=l 


(60) 


The  erased  RSO  (erasure  due  to  a  star  streak  mask)  will  not  invalidate  the  detection.  However, 
when  this  situation  occurs,  it  must  be  reported  and  the  concerned  frames  (where  star  streak  or 
streaks  overlap  the  RSO)  must  not  be  considered  for  the  RSO  brightness  estimation.  Concretely, 
if  there  are  9  images  used  for  RSO  detection,  and  in  two  of  the  nine  images  star  streaks  overlap 
the  RSO,  only  the  seven  other  images  (non  star-streak/RSO  overlap  images)  are  used  in  the  RSO 
brightness  calculations. 
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Table  I: 

Result  of  the  tests  for  an  RSO  with  a  SNR=4  and  a  sequence  of  9  images 

fS('J)  >9a„ 

Does  the  global  test 
succeed? 

I  °(ij)  >2an 

Do  the  individual  tests 
succeed? 

rr  ( i_D(io)  - 2a") 

k= 1 

Do  the  K  individual 
tests  succeed? 

RSO  +  noise 

SNR=4  in  individual 
frame,  SNR=12  with 
the  summation  of  9 
frames. 

Always 

Always 

Always 

Noise  only 

Never 

Sometimes; 

noise  events  above 

2  an  (3%  of  pixels). 

Never 

Cosmic-ray  hit  + 
noise 

Often; 

bright  flash 

Sometimes: 

the  frame  that 
contains  the  cosmic- 
ray  hit  and  noise 
events  above  2an. 

Never 

Noise  with 
background  residue 

Never 

Sometimes; 

noise  events  + 
background  residue 
above  2 cn 

Possible 

5.8  Detection  of  a  RSO  with  a  corruption  by  a  CR 

A  CR  will  not  prevent  the  detection  of  an  RSO.  However,  it  will  corrupt  the  calculation  of  the 
RSO  centroid  in  the  frame  where  it  appears.  Thus,  the  Cosmic  Ray  noise  artefact  must  be 
detected  and  the  corrupted  centroid  must  be  ignored.  When  an  RSO  is  detected  with  the  previous 
tests,  the  intensity  regularity  can  be  verified.  For  a  detection  ‘Ay’  at  coordinates  (i,j),  the  expected 
RSO  intensity  ‘<iRS0>’  js  Well  estimated  with  the  median  filter  (applied  on  the  drift  compensated 
images): 

For  Ay=l:  <IRS0(i,j)>  =  FDmed(i,j)  =  median  (I  !,,(i,j).  FD2(i,j)...  ,  FDK(i,j)).  (61) 

The  accepted  brightness  increment,  for  a  specific  observation  ‘k’,  should  be  less  than  a  predefined 
tolerance  ‘ht’; 

I  Dk(iJ)  -  <IRS0(i,j)>  <  ht .  (62) 

The  variation  tolerance  should  be  something  close  to  2an  (maybe  3  an)  for  faint  RSO  or  a  fraction 
of  the  brightness  for  brighter  RSO.  Only  the  photon  and  sensor  noises  should  be  responsible  for 
the  variation,  perhaps  including  the  PSF  influence  and  possibly  a  certain  tolerance  to  the  jitter. 
‘ht’  could  be  defined  with  an  ad-hoc  function  similar  to: 
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ht  =  3  sqrt  (an2  +  c  <IRS0(i,j)>  ), 


(63) 


which  considers  the  sensor  noise  and  Poisson  nature  of  the  signal.  This  assures  a  minimum 
tolerance  that  depends  on  the  noise  level.  This  tolerance  increases  with  stronger  RSO  signals 
(where  variations  will  be  greater).  The  function  takes  into  account  the  Poisson-distribution  nature 
of  the  photonic  signal.  The  constant  ‘c’  represents  the  ratio  of  photon-to-DN  (digital  number) 
conversion.  One  kind  of  error  that  Eq.  63  does  not  take  into  account  is  the  signal  quantization  or 
‘pixelization’. 

If  one  of  the  RSO  measurements  is  much  brighter  than  the  predefined  tolerance  for  an  individual 
image,  then  a  cosmic-ray  hit  can  be  suspected.  Recall  that  the  star  streaks  have  already  been 
erased  and  that  the  test  of  Eq.  62  does  not  take  intensity  decrement  into  account.  Thus,  when  the 
product  of  all  intensity  tests  succeeds  and  the  brightness  is  very  constant,  then  an  uncorrupted 
RSO  detection  can  be  declared  with; 

K 

1  =  ]~[  [  {  |  TDk(i,j)  -  <IRS0(i,j)>}  <  h,  ]  ;  constant  brightness  (64) 

k= 1 

Else,  when  one  of  the  ‘k’  tests  fails,  it  multiplies  the  result  by  zero.  This  means  that  a  pixel  is 
much  brighter  than  the  tolerated  variation.  This  is  an  indication  of  a  corruption  by  a  cosmic -ray 
hit.  In  that  case  the  measurement  must  be  ignored  in  the  evaluation  of  the  RSO  position. 

Note  that  this  test  does  not  work  with  a  fast  tumbling  satellite  -  that  implies  significant  brightness 
variation  over  time.  However,  the  tumbling  is  supposed  to  be  an  a  priori  knowledge  when  a 
satellite  acquisition  is  tasked.  Thus  when  such  a  case  is  present,  it  will  require  different  detection 
conditions  and  potentially  algorithmic  customizations  to  account  for  significant  variations  in  RSO 
signal  strength  during  the  RSO  imaging  sequence. 

5.9  Refinements  of  detection  algorithms  for  the  improvement 
of  performance  and  stability 

It  is  possible  to  improve  the  detection  performance  with  a  few  ad-hoc  modifications.  These 
modifications  take  into  account  that  1)  processing  artefacts  can  be  attenuated,  2)  the  signal  is  not 
completely  exploited  and  3)  evaluation  errors  can  be  corrected  or  compensated.  Here  are  a  few 
examples. 

5.9.1  Elimination  of  the  background  residue 

When  the  background  is  estimated  with  the  iterative  method  of  Eq.  36,  it  is  expected  that  it  leaves 
a  residue  \‘/2f  that  is  almost  null.  However,  when  the  image  sequence  is  summed  into  a  single 
image  containing  the  total  signal  (Eq.  53),  the  sum  of  the  residues  may  create  a  significant  signal. 
Hence,  Eq.  36  can  be  applied  again  on  the  image  of  total  signal  ‘I1’,  evaluating  the  global  residue, 
which  creates  a  new  background.  It  is  not  expected  to  be  more  than  a  fraction  of  the  noise  level 
multiplied  by  the  number  of  image  frames.  This  value  is  small  but  probably  measurable  and  it 
can  affect  the  detection  tests.  This  new  global  background  can  be  subtracted,  reducing  the  total 
residue  and  improving  the  detection  condition  of  Eq.  54. 

<Bf>  =  Is  0  Targe  impulse  response’.  (65) 
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Is'  =  Is  -  <BZ>  (66) 

A  more  sophisticate  method  could  consider  recalculating  the  image  summation  with  exclusion  of 
signal  above  the  noise  level. 

5.9.2  Integration  over  the  PSF  area 

The  previous  algorithms  do  not  take  into  account  that  the  signal  can  be  degraded  by  the  PSF. 
Detection  sensitivity  can  be  increased  by  integrating  the  signal  over  the  effective  PSF  area.  This 
can  be  achieved  simply  by  using  local  summation.  According  to  Annex  B  of  Ref.  2,  the  SNR  is 
optimal  when  the  summed  area  matches  the  PSF  area.  For  example,  when  the  PSF  width  is  3 
pixels  at  the  half  height,  a  summation  over  an  area  of  3x3  pixels  can  be  used.  This  can  effectively 
double  the  SNR.  Flere,  the  ‘square  root’  rule  cannot  be  used  because  the  pixels  do  not  all  have 
the  same  weight.  This  number  (double  SNR  for  a  3x3  window)  was  obtained  by  integrating  a 
Gaussian  PSF  profile  (3-pixel  width  at  half  height)  in  the  presence  of  Gaussian  noise. 

For  a  detection  with  integrated  area,  one  has  only  to  replace  the  images  ‘Is’  and  TDk’  in  Eq.  58  by 
their  equivalent  locally  summed  images  like  (here  for  a  PSF  width  equal  to  2n+l  pixels); 


i+n  j+n 

I-DkPSF(i,j)=  Z  Z  I  Dk(i,j) 

i—n  j—n 


i+n  j+n 

IZPSF(iJ)=  Z  Z  I  (i.j) 

i—n  j—n 


(67  A) 

(67B) 


5.9.3  Increase  the  tolerance  to  drift  estimation  errors 

The  detection  algorithms  also  consider  that  the  images  are  perfectly  realigned.  This  may  not  be 
the  case.  The  astrometric  calibration  of  the  celestial  pointing  vector  ‘PA’  may  leave  an  error  of 
one  or  two  pixels  in  the  drift  estimation.  The  RSO  may  not  be  perfectly  realigned  in  the  corrected 
image  sequence. 

The  detection  condition  of  Eq.  58  can  also  be  slightly  modified  to  accommodate  a  certain  level  of 
error  in  the  image  alignment.  Rather  than  comparing  aligned  pixels,  it  may  be  better  to  compare 
local  maxima.  The  detection  algorithm  could  accept  a  situation  like  this  one:  a  target  succeeds 
the  tests  in  the  successive  frames  at  the  locations:  (i,j)i,  (i.j+l)2,  (i-1 0)3,  etc.,  with  more  or  less 
one  pixel  of  jitter  from  frame  to  frame.  With  appropriate  threshold  values  (and  background 
correctly  removed)  there  is  almost  no  chance  that  noise  pattern  succeeds  in  generating  a  false 
alarm,  even  though  this  probability  increases.  A  simple  method  can  be  implemented  to  increase 
position  error  tolerance;  simply  replace  images  TDk’  by  one  of  local  maxima  such  as; 

TDkmax(i,j)  =  max  (I'Dk(i-n:i+n,j-n:j+n),  (68) 

where  ‘n’  is  the  tolerated  displacement  and  ‘i-n:i+n’  the  index  sequence 
‘i-n,  ...,  ,  i,  ...,  i+n’  (like  the  Matlab  notation).  This  algorithm  modification  can  be  chained  in 
cascade  with  Eq.  67. 
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5.10  Calculation  of  RSO  centroid 


Once  the  detections  are  known,  the  exact  position  of  the  RSO  must  be  determined.  This 
calculation  is  done  with  all  the  pixels  that  may  contain  signal,  even  those  that  do  not  contribute  to 
the  detection  but  surround  the  detected  pixels.  Thus,  a  small  sub-image 
TV,  centered  on  the  detected  pixel  ‘Ay’  and  larger  than  the  PSF,  is  extracted; 


I"Dkn  =  extract_kemel  (I'Dk,  Al|  n)  (69) 

The  ‘extractkemel’  is  a  function  that  selects  a  group  of  pixels  centered  on  the  nth  detection  at 
location  (i,j)  with  a  width  determined  by  the  PSF  width.  The  center  of  intensity  (same  calculation 
as  the  center  of  mass)  is  calculated  with  a  sub-pixel  accuracy  and  this  location  is  reported  in  the 
frame  of  reference  of  the  image  TDk’.  Before  the  calculation  of  sidereal  coordinates  the  drift 
value  ‘Xk’,  which  was  removed  by  the  image  realignment,  must  be  re-introduced  into  the  final 
position; 

RAkn  =  (i,j)RA  =  Xk  +  center  of  intensity  (l"Dkn).  (70) 

So  the  astrometric  conversion,  known  for  this  frame,  can  be  applied. 

RknA(tk)  =  @k-'  (  n  -'(RAkn))  (71) 

5.11  Cleaning  flashes  caused  by  cosmic-ray  hits 

In  some  cases,  it  might  be  interesting  to  remove  the  cosmic-ray  hits  from  the  images,  even  though 
this  was  not  necessary  for  the  detection.  In  fact,  it  is  usually  believed  that  the  cosmic-ray  hits 
must  be  erased  first.  The  acquisition  of  an  image  sequence  is  normally  done  for  this  purpose  and 
a  ‘voting’  filter  eliminates  the  bright  pixels  that  appear  in  only  one  frame.  If  this  is  done  first,  this 
would  also  erase  the  star  streaks  (which  appears  at  different  places  in  time)  and  make  astrometry 
calibration  difficult.  This  would  also  erase  the  RSO  if  there  is  pointing  drift.  So,  it  is  better  to 
leave  them.  But  they  can  be  erased  in  the  final  image  for  cosmetic  or  human  viewing  purposes. 
A  filter,  similar  to  Eq.  60,  can  be  used  to  declare  explicitly  the  cosmic-ray  hits  when  a  pixel  is 
brighter  than  the  expected  value  (plus  the  tolerance  variation).  Eq.  62  can  be  rewritten  for  CR 
detection,  i.e.; 

If  [  (I“Dk(i,j)  -  Imed(i,j)  )  >  ht  ]  ;  this  is  a  cosmic-ray.  (72) 

In  such  a  case,  the  cosmic  ray  can  be  erased  by  replacing  the  pixel  value  by  the  most  probable 
value  given  by  the  median  filter; 

FDCk(ij)  =  Imed(i,j)  (73) 

Td"c’  is  read;  the  drift-compensated  0’°’)  and  cosmic-ray-hit  cleaned  C'c’)  image. 
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5.12  Detecting  serendipitous  RSOs. 

The  previous  processing  framework  is  designed  for  the  detection  of  expected  RSOs  with  the 
maximum  possible  sensitivity.  The  drawback  of  the  method  is  that  it  does  not  have  the  capability 
of  detecting  unexpected  RSOs,  even  bright  ones.  This  is  a  shortfall  of  methods  designed  with  a- 
priori  knowledge  such  as  matched  fdter  and  TRM  acquisition. 

For  detection  of  serendipitous  RSOs  acquired  in  SSM  (sidereal  stare  mode),  it  is  recommended  to 
use  other  detection  techniques  such  as  those  described  in  Refs.  7  and  8.  Basically,  these  methods 
detect  candidates  in  individual  frames  and  then  look  for  objects  that  have  coherent  alignment  in 
the  image  sequence.  Obviously,  only  bright  objects  can  be  detected  in  each  frame  of  the 
sequence  and  be  noted.  These  methods  are  not  as  sensitive  as  the  method  presented  in  this 
document. 

In  TRM,  when  an  unexpected  RSO  is  present  in  the  FOV,  there  are  two  possible  cases.  It  may 
have  a  motion  very  close  to  the  other  expected  object  and  it  can  be  detected  as  a  member  of  the 
group  of  RSOs  (see  Section  4.9).  The  other  possibility  is  that  it  has  a  very  different  motion  and 
appears  like  a  streak  in  images.  If  the  streak  alignment  is  the  same  as  the  stars,  it  will  be  detected 
as  a  different  star  in  each  frame.  The  only  way  to  detect  it  is  to  observe  that  there  is  an 
unidentified  star  in  each  frame  and  the  position  pattern  is  coherent  in  the  image  sequence.  This 
would  require  ad-hoc  algorithms  which  are  not  yet  developed.  If  the  alignment  is  different,  it 
will  be  equivalent  to  a  cosmic-ray  hit  (CR)  in  the  individual  frames.  However,  it  creates  a 
coherent  aligned  pattern  in  the  non-realigned  image  sequence.  Here  is  an  algorithm  framework 
that  could  detect  bright  serendipitous  RSOs: 

1  -  detect  and  remove  stars  in  each  frame, 

2-  sum  all  the  image  frames  into  a  single  global  image, 

3-  detect  streak  patterns  (whose  orientation  is  different  than  the  stars)  with  Hough  transform, 

4-  look  for  chopped  streaks. 

This  will  detect  all  CRs  and  all  bright  serendipitous  RSOs.  The  difference  between  these  two 
groups  of  object  is  that  CRs  are  single  uniform  streaks  (with  various  lengths)  while  RSOs  will 
appear  as  chopped  streaks  (several  exposures  with  dead  time  between  the  exposures).  The  regular 
chopping  can  be  used  for  RSO  detection.  This  may  be  an  interesting  algorithmic  framework  to 
further  develop. 
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6  Processing  scheme  for  RSO  detection 


For  the  detection  of  satellites  in  TRM  mode,  the  global  processing  scheme  is  presented  in  the 
following  figures.  This  scheme  is  based  on  the  image  model  developed  in  the  previous  chapters. 
It  illustrates  how  algorithms  are  sequenced  along  with  their  dependencies.  For  detection  in 
sidereal  tracking  mode,  the  processing  methodology  has  already  been  published  in  Refs.  1  and  2. 

In  this  processing,  the  image  background  is  the  first  thing  to  be  eliminated.  Section  5.2  explained 
that  it  can  be  estimated  by  different  methods.  Once  it  is  adequately  estimated,  it  is  subtracted 
from  the  image,  producing  a  background-free  image  T8’. 

The  next  step  consists  of  detecting  the  stars.  This  uses  the  satellite  pointing  and  attitude 
parameters  (pointing  coordinates,  motion  and  camera  rotation)  for  the  estimation  of  the  length 
and  orientation  of  the  star  streaks  in  the  image.  This  will  produce  the  convolution  kernel 
‘Z(PA(t)|t,t+dt)’  used  by  the  iterative  matched  filter.  If  there  are  enough  very  bright  stars,  this 
method  can  be  replaced  by  a  faster  algorithm  (a  basic  threshold  method  and  Hough  transforms) 
but  the  matched  filter  method  guarantees  a  maximum  of  detected  stars.  Once  the  stars  are 
detected  (‘I+s’),  they  are  extracted  and  their  centres  of  intensity  are  evaluated.  Then,  the  list  of 
star  relative  brightness  and  location  ‘Sk(n,i,j,an)’  is  built. 

Before  performing  the  RSO  detection,  it  is  preferable  to  remove  the  star  streaks  from  the  images. 
This  makes  the  evaluation  of  the  total  RSO  signal  more  accurate  and  also  decreases  the  number  of 
false  alarms.  The  bright  star  streaks  (and  possibly  every  streak)  are  already  detected  (‘I+s’). 
Their  binary  extraction  masks  ‘Msk’  are  also  already  created.  They  can  be  used  to  reset  the  streak 
to  zero  in  the  images  (‘Fs’).  Because  of  the  detection  method  that  compares  every  frame  of  the 
sequence  for  the  detection  of  the  static  object,  it  is  not  certain  that  the  star  erasing  is  really 
required.  Tests  will  have  to  be  performed  to  determine  if  streak  removal  improves  detection 
performance.  On  the  contrary,  if  tests  demonstrate  that  the  star  streak  must  be  erased  and  that  the 
‘reset’  method  is  not  accurate  enough,  Ref.  1  explains  a  more  accurate  method  that  removes  stars 
without  leaving  artifacts.  This  method  could  be  adapted  for  the  current  case  but  this  is  probably 
not  necessary. 

Once  the  stars  are  detected,  the  images  can  be  calibrated  astrometrically,  i.e.,  by  establishing  the 
relation  between  the  pixel  coordinates  and  their  equivalent  sidereal  (Ra,  Dec)  coordinates.  The 
star  pixel  locations  are  first  corrected  for  the  geometric  distortion  caused  by  the  optics.  Then,  this 
information  is  compared  with  a  star  database,  using  commercial  software  (like  ‘PinPoint  or 
‘TheSky’).  This  determines  the  relation  ‘@k’  between  the  pixel  indexes  and  sidereal  coordinates. 
This  calibration  also  determines  the  exact  telescope  pointing  direction  ray  ‘PA(t)\  This  also 
facilitates  the  evaluation  of  pointing  errors  and  drifts. 

Once  the  drift  is  known,  the  image  frames  Ts’  are  re-aligned  TD’  and  the  RSO  appears  at  exactly 
the  same  place  in  every  frame.  This  allows  the  summation  of  all  images  ‘I1’,  which  improves  the 
SNR  and  increases  global  detection  sensitivity.  Sensitivity  can  also  be  increased  by  integrating 
the  RSO  signal  over  the  PSF  area  ‘IIPSF’  and/or  using  the  nearest  maxima  ‘pDkmax,_  Finally,  the 
RSO  detection  ‘Ay’  is  performed  with  the  double  threshold  test:  1)  the  global  detection  with  the 
high  threshold  value  and  2)  the  product  of  all  individual  detections  with  the  lower  threshold 
value. 
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Where  an  RSO  is  detected,  the  brightness  is  evaluated  and  the  brightness  constancy  is  verified.  If 
the  brightness  suddenly  increases  for  an  individual  frame,  there  is  a  high  probability  that  this 
detection  is  corrupted  by  a  CR  or  a  star  streak.  The  detection  with  the  CR  contamination  can  be 
ignored  in  the  detection  report,  or,  at  least,  in  the  brightness  estimation.  The  star  streak  proximity 
can  also  be  checked  with  the  realigned  streak  mask  ‘MS  Dk’. 

When  detection  occurs,  it  must  be  reported.  However,  the  pixel  indices  need  to  be  rewritten  in 
the  original  frame  coordinates  by  cancelling  the  image  realignment  offset.  The  cumulative  drift, 
which  was  subtracted  from  the  image  coordinates,  must  now  be  reintroduced  in  the  detection 
coordinates.  Afterward  the  geometric  correction  ‘ft"1’  must  be  applied  to  these  coordinates. 
Finally,  it  remains  only  to  apply  the  astrometric  conversion  ‘@k-1’  to  obtain  the  real  RSO 
coordinates. 
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Figure  4:  Processing  scheme  for  the  detection  ofRSO  acquire  in  TRM  mode 
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Figure  5:  Processing  scheme  for  the  detection  of  RSO  acquire  in  TRM  mode: 

continued from  Fig  4 
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Figure  6:  Processing  scheme  for  the  detection  of  RSO  acquire  in  TRM  mode: 

continued from  Fig  5 
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7  Conclusion 


The  strong  expertise  acquired  with  the  past  development  of  RSO  detection  algorithms  (Refs.  1  to 
4)  allows  the  design  of  these  algorithms  which  are  specifically  adapted  for  the  detection  in 
satellite  tracking  mode.  These  algorithms  have  not  yet  been  fully  coded  nor  tested,  but  the  most 
challenging  components  have  already  been  tested  in  previous  works  (Refs.  1  and  2).  They  are  the 
iterative  matched  filter,  the  background  and  star  removal  methods,  the  astrometric  calibration, 
image  masking  and  segmentation,  etc.  All  these  algorithms  are  recycled  and  rearranged  for  the 
current  case.  This  new  method  takes  advantage  of  an  acquisition  method  that  acquires  an  image 
sequence,  rather  than  a  single  one,  which  makes  detection  easier,  more  sensitive  and  more 
reliable.  This  introduces  the  need  to  measure  the  drift  and  realign  the  image,  but  this  does  not 
represent  a  difficult  technical  challenge. 

This  document  will  be  used  to  guide  the  software  development  that  will  be  required  by  space- 
based  systems  under  development  like  NEOSSat  or  SAPPHIRE.  Some  algorithms  may  be 
changed  for  various  reasons  (performance,  accuracy,  fewer  processing  artifacts,  sensitivity,  etc), 
but  the  global  approach  should  remain  similar.  Tests  will  have  to  be  performed  to  fine  tune  the 
threshold  values,  for  the  evaluation  of  the  detection  sensitivity  and  for  the  drawing  of  the 
probability  of  detection  graphs.  All  these  tests  will  have  to  be  performed  with  simulated  data  first 
because  the  real  data  may  not  be  available  at  that  time.  The  MOST  satellite  could  provide  similar 
real  data,  but  this  satellite  was  not  conceived  with  a  similar  capability  (like  NEOSSat)  for  the 
acquisition  of  fast  image  sequences  (short  exposure  and  short  download  delay).  Possibly,  images 
from  the  MSX/SBV  sensor  (Space  Based  Visual)  or  the  new  SBSS  (Space  Based  Surveillance  of 
Space)  satellite  could  be  used  for  testing  with  real  data  if  this  data  is  made  available  to  the 
developer  team. 

These  algorithms  were  developed  in  the  context  of  a  space-based  surveillance  system  but  they  can 
also  be  used  by  ground  stations  when  they  acquire  images  in  TRM  mode.  The  only  characteristic 
of  a  space-based  sensor  is  that  it  is  exposed  to  cosmic  rays  and  the  images  suffer  from  a  source  of 
noise  that  is  not  ordinarily  present  in  a  ground  station  data.  On  the  other  hand,  data  are  not 
affected  by  the  atmospheric  transmission  interferences,  including  turbulence.  In  TRM  mode,  both 
pointing  systems  are  moving  and  only  the  final  apparent  RSO  motion  is  taken  into  account;  the 
fact  that  one  sensor  is  moving  while  the  other  is  at  a  fixed  location  is  not  an  important  parameter 
for  image  processing  algorithms. 

This  method  should  be  very  sensitive.  It  is  expected  to  be  able  to  detect  RSO  as  faint  as  SNR>3. 
It  uses  the  fact  that  the  fusion  of  the  image  sequence  provides  an  improved  SNR.  At  the  same 
time,  single  detection,  in  the  individual  frames,  can  tolerate  a  very  high  rate  of  false  alarms.  The 
probability  that  there  are  false  alarms  at  the  same  coordinates  in  every  frame  is  almost  null. 

Contrary  to  popular  opinion,  cosmic-ray  hits  (CR)  do  not  have  to  be  erased  in  the  images. 
Usually,  the  first  processing  step  is  a  voting  system  in  whuch  the  flashes  cause  by  CRs  are 
detected  and  erased.  In  TRM  mode,  this  technique  would  affect  the  star  streaks,  which  have  the 
same  behavior  as  the  CRs.  With  current  detection  algorithms,  CRs  are  ignored;  when  a  CR  trigs 
an  alarm,  it  is  detected  in  a  single  frame,  not  the  whole  sequence.  In  addition,  a  real  detection  can 
be  corrupted,  in  one  or  several  frames,  by  CRs  or  star  streaks.  However,  when  a  possible 
detection  is  notified,  brightness  regularity  is  verified  and  the  frames  whose  brightness  suddenly 
becomes  excessively  hight  are  ignored.  The  keys  for  such  detection  performance  are:  the  number 
of  images  in  the  sequence  and  the  ability  to  realign  the  image  which  compensates  for  the  RSO 
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drift.  This  drift  compensation  is  possible  because  of  the  astrometric  calibration  and  the 
knowledge  about  the  prediction  of  the  RSO  position. 

This  processing  scheme  is  not  intended  for  on-board  real-time  processing;  it  requires  too  much 
processing  capability.  It  is  designed  for  post-acquisition  processing  by  ground  stations  where  all 
the  required  computing  resources  are  available.  Thus,  in  the  intended  mode  of  operations, 
imaging  data  are  downloaded  and  updated  RSO  orbital  parameters  calculated  on  the  ground.  This 
greatly  reduces  risks  for  algorithm  implementation  as  revisions  (if  necessary)  can  be  readily 
implemented  (no  remote  software  updates  required).  Also,  it  permits  on-ground  diagnostics  for 
cases  when  an  RSO  is  not  detected  or  ambiguous  detection  results  occur. 
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List  of  symbols/abbreviations/acronyms/initialisms 


CR 

Cosmic-ray  hit. 

RSO: 

Space  Resident  Object  (satellite). 

TLE: 

‘Two-Lines  Element  set’:  SSN’s  format  for  the  parameters 
used  to  calculate  the  orbital  position  of  a  satellite. 

SSN: 

Space  Surveillance  Network. 

FOV: 

Field  of  view. 

PSF: 

Point  Spread  Function,  which  is  also  the  optical  impulse 
response. 

General: 

T: 

Temperature. 

t: 

Time. 

t  (or  dt): 

Exposition  time. 

® 

Convolution  operator. 

Rot(t): 

Rotation  of  the  sensor  at  time  ‘t\ 

Coordinates  systems: 

(ij)k: 

pixel  index  (i,j)  of  image  frame  ‘k\ 

(Ra,Dec,t): 

(Right  ascension.  Declination,  time). 

k  (lower  case): 

Frame  index 

K  (upper  case): 

Number  of  frames  in  the  image  sequence. 

Coordinate  transformations: 

@k- 

Coordinate  transformation;  (i”,j”)  =  @k(Ra,Dec) 

ft: 

Geometric  distortion  (optical  distortion);  IAk  (i',j’)  =  l\  ft  (i 

0: 

Optical  PSF;  IAk(i,j)  =  0  (IAk  (i’j’)) 

Scene  elements: 

Astronomical  scene. 

an. 

Intensity  of  the  nth  star. 

(Ra,Dec)n: 

aR:  Intensity  of  the  RSO. 

Location  of  the  nth  star. 

(Ra,Dec)RA: 

Location  of  the  RSO. 
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Z(PA(t)|t,t+dt): 

Streak  shape  produced  by  the  pointing  motion  PA(t)  during  the 
time  interval  t,t+dt. 

Z((RA(t)-PA(t))|t,t+dt) 

Streak  shape  produced  by  the  drift. 

Position  and  pointing: 

Rp(TLE,t)  =  (Ra,  Dec,  t)Rp: 

Predicted  position  of  the  satellite. 

RA(t)  =  (Ra,  Dec,  t)RA  : 

Real  position  of  the  satellite 

P°(t)  =  (Ra,  Dec,  t)Rp: 

Pointing  vector,  determined  with  the  onboard  instruments. 

PA(t)  =  (Ra,  Dec,  t)PA: 

Real  pointing  vector,  determined  with  astrometry. 

D(t)  : 

Drift  between  two  consecutive  image  frames. 

Fp(t)  : 

Pointing  offset. 

FR(t)  : 

RSO  position  offset. 

Images  components: 

IA: 

Image  component  containing  only  the  stars. 

IR: 

Image  component  containing  only  the  RSO. 

1°: 

Original  (acquired)  image. 

IB: 

Background-free  image. 

I+s: 

Image  component  containing  only  the  star  streaks. 

I+Skn: 

A  sub-image  containing  only  the  nth  streak  of  the 
kth  image  frame. 

Is: 

Image  frame  where  star  streaks  (and  also  the 
background)  are  removed. 

I-D: 

Image  frame  where  the  drift  is  compensated 
(image  realigned). 

f: 

Summation  of  the  realigned  image  sequence. 

B: 

Image  background. 

<B>: 

Estimated  image  background. 

SI: 

Background  residue  (after  background  subtraction). 

n: 

noise 

cn: 

noise  standard  deviation. 

Ck(0  : 

Image  component  containing  the  cosmic-ray  hits  acquired 
during  the  exposition  time. 

-A*: 

Convolved  image;  the  convolution  peaks  indicate  the 
streak  locations. 

5): 

Dark  frame. 
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Lj(PA(tk),  t): 

L2(PA(tk),  t): 
Ms: 

Sk.„: 

Sk(n,i,j): 


Astronomical  scene  background  (zodiacal  light, 
diffuse  nebulae,  etc.). 

Image  background  component  caused  by  stray  light. 
Binary  segmented  mask  image  used  to  extract  streaks. 
Sub-image  containing  only  one  streak. 

List  of  the  ‘n’  star  centroid  coordinates  (  i,j). 


Detection: 

hh: 

hi: 


ht: 


High  level  detection  threshold;  for  global  detection. 

Low  level  detection  threshold;  for  detection  in  individual 
frame. 

Brightness  tolerance  threshold;  for  cosmic-ray  hit 
detection. 


Aijk 

<IRso> 


Detection  at  coordinates  (i,j)  in  the  frame  ‘k’. 
Once  detected,  estimation  of  the  RSO  brightness. 
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In  the  context  of  the  surveillance  of  space,  known  resident  space  objects  (RSO),  i.e.  active 
satellites  and  space  debris,  need  to  be  acquired  periodically  to  ensure  that  the  knowledge  of  their 
orbital  parameters  is  up-to-date.  The  most  sensitive  observation  method  consists  of  acquiring  an 
image  sequence  while  the  target  is  actively  tracked  by  the  optical  system  such  the  RSO  appears 
stationary  in  the  image  sequence.  This  generates  an  incredible  amount  of  data  which  requires 
efficient  detection  algorithms  and  automated  processing  software.  In  this  context,  this  document 
presents  models  for  image  formation,  acquisition,  and  processing.  From  the  observation  request 
to  the  production  of  the  images,  the  stream  of  events  is  briefly  modeled.  This  ensures  that  a- 
priori  knowledge  and  strategy  will  be  applied  to  efficiently  process  the  acquired  data.  These 
models  produce  a  series  of  algorithms  used  to  separate  the  image  components  and  detect  the 
RSO.  This  document  presents  algorithm  diagrams  that  can  be  used  to  develop  the  processing 
software. 

Dans  le  cadre  de  la  surveillance  de  I’espace,  les  objets  spatiaux  connus  en  orbite  (OSO),  i.e., 
satellites  actifs  ou  debris  spatiaux,  ont  besoin  d’etre  acquis  periodiquement  pour  s’assurer  que  la 
connaissance  de  leurs  parametres  orbitaux  soit  maintenu  a  jour.  La  methode  d’observation  la 
plus  sensible  consiste  a  acquerir  une  sequence  d’images  lorsque  la  cible  est  poursuivie 
activement  par  le  systeme  optique  et  que  I’OSO  apparaTt  stationnaire  dans  la  sequence 
d’images.  Cela  genere  une  quantite  incroyable  de  donnees  qui  requierent  des  algorithmes  de 
detection  efficaces  et  des  logiciels  de  traitement  automatise.  Dans  ce  contexte,  ce  document 
presente  des  modeles  pour  la  formation,  I’acquisition  et  le  traitement  des  images.  De  la  requete 
d’observation  a  la  production  des  images,  la  chaine  des  evenements  est  modelisee  brievement. 
Cela  permet  de  voir  les  connaissances  et  strategies  a  priori  qui  peuvent  etre  utilisees  pour  traiter 
efficacement  le  donnees  acquises.  Ces  modeles  produisent  une  serie  d’algorithmes  utilises  pour 
separer  les  composantes  de  I’image  et  detecter  I’OSO.  Ce  document  presente  les  diagrammes 
d’algorithmes  qui  seront  utilises  pour  developper  le  logiciel  de  traitement. 
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